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Abstract. A preconditioning theory is presented which establishes sufficient condi- 
tions for multiplicative and additive Schwarz algorithms to yield self-adjoint positive 
definite preconditioners. It allows for the analysis and use of non-variational and non- 
convergent linear methods as preconditioners for conjugate gradient methods, and it is 
applied to domain decomposition and multigrid. It is illustrated why symmetrizing may 
be a bad idea for linear methods. It is conjectured that enforcing minimal symmetry 
achieves the best results when combined with conjugate gradient acceleration. Also, it is 
shown that absence of symmetry in the linear preconditioner is advantageous when the 
linear method is accelerated by using the Bi-CGstab method. Numerical examples are 
presented for two test problems which illustrate the theory and conjectures. 



Domain decomposition (DD) and multigrid (MG) methods have been studied exten- 
sively in recent years, both from a theoretical and numerical point of view. DD methods 
were first proposed in 1869 by H. A. Schwarz as a theoretical tool in the study of ellip- 
tic problems on non-rectangular domains [22] . More recently, DD methods have been 
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reexamined for use as practical computational tools in the (parallel) solution of general 
elliptic equations on complex domains [16]. MG methods were discovered much more 
recently [10]. They have been extensively developed both theoretically and practically 
since the late seventies [6, 11], and they have proven to be extremely efficient for solving 
very broad classes of partial differential equations. Recent insights in the product nature 
of certain MG methods have led to a unified theory of MG and DD methods, collectively 
referred to as Schwarz methods [5, 9, 27]. 

In this paper, we consider additive and multiplicative Schwarz methods and their accel- 
eration with Krylov methods, for the numerical solution of self-adjoint positive definite 
(SPD) operator equations arising from the discretization of elliptic partial differential 
equations. The standard theory of conjugate gradient acceleration of linear methods re- 
quires that a certain operator associated with the linear method - the preconditioner - 
be symmetric and positive definite. Often, however, as in the case of Schwarz-based 
preconditioners, the preconditioner is known only implicitly, and symmetry and positive 
definiteness are not easily verified. Here, we try to construct natural sets of sufficient 
conditions that are easily verified and do not require the explicit formulation of the pre- 
conditioner. More precisely, we derive conditions for the constituent components of MG 
and DD algorithms (smoother, subdomain solver, transfer operators, etc.), that guarantee 
symmetry and positive definiteness of the preconditioning operator which is (explicitly 
or implicitly) defined by the resulting Schwarz method. 

We examine the implications of these conditions for various formulations of the stan- 
dard DD and MG algorithms. The theory we develop helps to explain the often observed 
behavior of a poor or even divergent MG or DD method which becomes an excellent 
preconditioner when accelerated by a conjugate gradient method. We also investigate 
the role of symmetry in linear methods and preconditioners. Both analysis and numer- 
ical evidence suggest that linear methods should not be symmetrized when used alone, 
and only minimally symmetrized when accelerated by conjugate gradients, in order to 
achieve the best possible convergence results. In fact, the best results are often obtained 
when a very nonsymmetric linear iteration is used in combination with a nonsymmetric 
system solver such as Bi-CGstab, even though the original problem is SPD. 

The outline of the paper is as follows. We begin in §2 by reviewing basic linear meth- 
ods for SPD linear operator equations, and examine Krylov acceleration strategies. In 
§3 and §4, we analyze multiplicative and additive Schwarz preconditioners. We de- 
velop a theory that establishes sufficient conditions for the multiplicative and additive 
algorithms to yield SPD preconditioners. This theory is used to establish sufficient con- 
ditions for multiplicative and additive DD and MG methods, and allows for analysis of 
non-variational and even non-convergent linear methods as preconditioners. A simple 
lemma, given in §5, illustrates why symmetrizing may be a bad idea for linear methods. 
In §6, results of numerical experiments obtained with finite-element-based DD and MG 
methods applied to some non-trivial test problems are reported. 

2. Krylov acceleration of linear iterative methods 

In this section, we review some background material on self-adjoint linear opera- 
tors, linear methods, and conjugate gradient acceleration. More thorough reviews can 
be found in [12, 18]. 

2.1. Background material and notation. Let H be a real finite-dimensional Hilbert 
space equipped with the inner-product (•, •) inducing the norm || • || = (•, -) 1 / 2 . H can be 
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thought of as, for example, the Euclidean space W n , or as an appropriate finite element 
space. 

The adjoint of a linear operator A E L("H,"H) with respect to (•, •) is the unique 
operator A T satisfying (Au, v) = (u, A T v) , Vu, v E H. An operator A is called self- 
adjoint or symmetric if A = A T ; a self-adjoint operator A is called positive definite or 
simply positive, if (Au,u) > , Vu E H, u ^ 0. If A is self-adjoint positive definite 
(SPD) with respect to (•, •), then the bilinear form (Au, v) defines another inner-product 

1 /2 

on %, which we denote as (-, -)a- It induces the norm || • ||a = (■, -)a ■ 

The adjoint of an operator M E L(H, H) with respect to (•, -)a, the A-adjoint, is the 
unique operator M* satisfying (Mu, v)a = (u, M*v)a , Vu, v E %. From this definition 
it follows that 

M* — A~ 1 M T A . (2.1) 

M is called A- self- adjoint if M = M* , and A-positive if (Mu, u) A > , Vit E H,u ^ 0. 

If N E L(Hi, H 2 ), then the adjoint of N, denoted as N T E L('H 2 , Hi), is defined as 
the unique operator relating the inner-products in Hi and H2 as follows: 

(Nu,v) n2 = (u,N T v) ni , VuEHi, VvEH 2 . (2.2) 

Since it is usually clear from the arguments which inner-product is involved, we shall 
often drop the subscripts on inner-products (and norms) throughout the paper, except 
when necessary to avoid confusion. 

We denote the spectrum of an operator M as a(M). The spectral theory for self- 
adjoint linear operators states that the eigenvalues of the self-adjoint operator M are real 
and lie in the closed interval [A m j n (M), A maa; (M)] defined by the Rayleigh quotients: 

,, r s ■ (Mu,u) . (Mu,u) 
X min (M) = mm \ X max M = max \ (2.3) 

Similarly, if an operator M is A-self-adjoint, then its eigenvalues are real and lie in the 
interval defined by the Rayleigh quotients generated by the A-inner-product. A well- 
known property is that if M is self-adjoint, then the spectral radius of M, denoted as 
p(M), satisfies p(M) = \\M\\. This property can also be shown to hold in the A-norm 
for A-self-adjoint operators (or, more generally, for A-normal operators [1]). 

Lemma 2.1. If A is SPD and M is A-self-adjoint, then p(M) = ||M||^. 

2.2. Linear methods. Given the equation Au = /, where A E L(H,H) is SPD, con- 
sider the preconditioned equation BAu = Bf, with B E L("H,"H). The operator B, 
the preconditioner, is usually chosen so that a Krylov or Richardson method applied to 
the preconditioned system has some desired convergence properties. A simple linear 
iterative method employing the operator B takes the form 

u n+l =y n_ BAu n + B f = (j _ BA)u n + £^ 

where the convergence behavior of (2.4) is determined by the properties of the so-called 
error propagation operator, 

E = I-BA. (2.5) 

The spectral radius of the error propagator E is called the convergence factor for the 
linear method, whereas the norm is referred to as the contraction number. We recall two 
well-known lemmas; see for example [17] or [20]. 

Lemma 2.2. For arbitrary f and u°, the condition p(E) < 1 15 necessary and sufficient 
for convergence of the linear method (2.4). 
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Lemma 2.3. The condition \\E\\ < 1, or the condition \\E\\a < L is sufficient for 
convergence of the linear method (2.4). 

We now state a series of simple lemmas that we shall use repeatedly in the following 
sections. Their short proofs are added for the reader's convenience. 

Lemma 2.4. If A is SPD, then BA is A- self -adjoint if and only if B is self-adjoint. 

Proof. Note that: (ABAu,v) = (BAu,Av) = (Au, B T Av). The lemma follows since 
BA = B T A if and only if B = B T . □ 

Lemma 2.5. If A is SPD, then E is A-self-adjoint if and only if B is self-adjoint. 

Proof. Note that: (AEu, v) = (Au,v)-(ABAu,v) = (Au,v)-(Au, (BA)*v) = {Au, (I- 
(BA)*)v). Therefore, E* = E if and only if BA = (BA)*. By Lemma 2.4, this holds if 
and only if B is self-adjoint. □ 

Lemma 2.6. If A and B are SPD, then BA is A-SPD. 

Proof. By Lemma 2.4, BA is A-self-adjoint. Also, we have (ABAu, u) = (BAu, Au) = 
(B 1/2 Au, B 1/2 Au) > , Vw ^ 0. Hence, BA is A-positive, and the result follows. □ 

Lemma 2.7. If A is SPD and B is self-adjoint, then \\E\\a = p{E). 

Proof. By Lemma 2.5, E is A-self-adjoint. By Lemma 2.1 the result follows. □ 

Lemma 2.8. If E* is the A-adjoint of E, then \\E\\\ = \\EE*\\a- 

Proof. The proof follows that of a familiar result for the Euclidean 2-norm [12]. □ 

Lemma 2.9. If A and B are SPD, and E is A-non-negative, then \\E\\a < 1. 

Proof. By Lemma 2.5, E is A-self-adjoint. As E is A-non-negative, it holds (Eu, u)a > 
0, or (BAu,u)a < (u,u)a- By Lemma 2.6, BA is A-SPD, and we have that < 
(BAu, u) A < (u, u) A , Vm ^ 0, which, by (2.3), implies that < A t < 1, VA» G cr(BA). 
Thus, p(E) = 1 — minj \ < 1. Finally, by Lemma 2.7, we have \\E\\a = p(E). □ 

We will also have use for the following two simple lemmas. 
Lemma 2.10. If A is SPD and B is self-adjoint, and E is such that: 

-Ci(u,u) A < (Eu,u)a < C 2 (u,u)a, Vm G H, 
for d > and C 2 > 0, then p(E) = \\E\\ A < max{Ci, C 2 }. 

Proof. By Lemma 2.5, E is A-self-adjoint, and by (2.3) X min (E) and X max (E) are 
bounded by —C\ and C 2 , respectively. The result then follows by Lemma 2.7. □ 

Lemma 2.11. If A and B are SPD, then Lemma 2.10 holds for some C 2 < 1. 

Proof. By Lemma 2.6, BA is A-SPD, which implies that the eigenvalues of BA are real 
and positive. Hence, we must have that \i(E) = 1 — \i(BA) < 1, VI Since C 2 in 
Lemma 2.10 bounds the largest positive eigenvalue of E, we have that C 2 < 1. □ 



SCHWARZ METHODS: TO SYMMETRIZE OR NOT TO SYMMETRIZE 



5 



2.3. Krylov acceleration of SPD linear methods. The conjugate gradient method was 
developed by Hestenes and Stiefel [13] as a method for solving linear systems Au = f 
with SPD operators A. In order to improve convergence, it is common to precondition 
the linear system by an SPD preconditioning operator B m A^ 1 , in which case the gener- 
alized or preconditioned conjugate gradient method results ([8]). Our goal in this section 
is to briefly review some relationships between the contraction number of a basic linear 
preconditioner and that of the resulting preconditioned conjugate gradient algorithm. 
We start with the well-known conjugate gradient contraction bound ([12]): 

\\e i+1 h < 2 (l - —^ 7 =] l|e°|U = 2 \\e°\\ A . (2.6) 

\ 1 + y/K A (BA) J 

The ratio of extreme eigenvalues of BA appearing in the derivation of the bound gives 
rise to the generalized condition number k a (BA) appearing above. This ratio is often 
mistakenly called the (spectral) condition number k(BA); in fact, since BA is not self- 
adjoint, this ratio is not in general equal to the usual condition number (this point is 
discussed in great detail in [1]). However, the ratio does yield a condition number in the 
A-norm. The following lemma is a special case of Corollary 4.2 in [1]. 



Lemma 2.12. If A and B are SPD, then 



n A {BA) = \\BA\\ A \\(BA)-i\\ A = . (2.7) 

Remark 2.13. Often a linear method requires a parameter a in order to be convergent, 
leading to an error propagator of the form E = I — aBA. Equation (2.7) shows that 
the A-condition number does not depend on the particular choice of a. Hence, one can 
use the conjugate gradient method as an accelerator for the method without a parameter, 
avoiding the possibly costly estimation of a good a. 

The following result gives a bound on the condition number of the operator BA in 
terms of the extreme eigenvalues of the error propagator E = I — BA; such bounds are 
often used in the analysis of linear preconditioners (cf. Proposition 5.1 in [26]). We give 
a short proof of this result for completeness. 

Lemma 2.14. If A and B are SPD, and E is such that: 

— Ci(u,u) A < (Eu,u) A < C 2 (u,u) A , Vu G K, (2.8) 
for Ci > and C 2 > 0, then the above must hold with C 2 < 1, and it follows that: 

ka(BA) < i±|L. 

Proof. First, since A and B are SPD, by Lemma 2.11 we have that C 2 < 1. Since 
(Eu, u) A = (u, u) A — (BAu, u) A , it is clear that 

(l-C 2 )(u,u) A < (BAu,u) A < (l + Ci)(u,u) A , Vm e H. 

By Lemma 2.6, BA is A-SPD. Its eigenvalues are real and positive, and lie in the interval 
defined by the Rayleigh quotients generated by the A-inner-product. Hence, that interval 
is given by [(1 — C 2 ), (1 + C\)\, and by Lemma 2.12 the result follows. □ 

Remark 2.15. Even if a linear method is not convergent, it may still be a good precon- 
ditioner. If it is the case that C 2 « 1, and if C\ > 1 does not become too large, then 
k a (BA) will be small and the conjugate gradient method will converge rapidly, even 
though the linear method diverges. 
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If only a bound on the norm of the error propagator E = I — BA is available, then the 
following result can be used to bound the condition number of BA. This result is used 
for example in [27]. 

Corollary 2.16. If A and B are SPD, and \\ I - BA\\ A < 5 < 1, then 

ka(BA) < (2.9) 
l — o 

Proof. This follows immediately from Lemma 2.14 with 5 = max{Ci, C2}. □ 

The next result connects the contraction number of the preconditioner to the contrac- 
tion number of the preconditioned conjugate gradient method. It shows that the conjugate 
gradient method always accelerates a linear method (if the conditions of the lemma hold). 

Lemma 2.17. If A and B are SPD, and \\I — BA\\ A < 5 < 1, then 5 cg < 5. 

Proof. An abbreviated proof appears in [27], a more detailed proof in [14]. □ 

2.4. Krylov acceleration of nonsymmetric linear methods. The convergence theory 
of the conjugate gradient iteration requires that the preconditioned operator BA be A- 
self-adjoint (see [2] for more general conditions), which from Lemma 2.4 requires that 
B be self-adjoint. If a Schwarz method is employed which produces a nonsymmetric 
operator B, then although A is SPD, the theory of the previous section does not ap- 
ply, and a nonsymmetric solver such as conjugate gradients on the normal equations [2], 
GMRES [21], CGS [23], or Bi-CGstab [25] must be used for the now non-A-SPD pre- 
conditioned system, BAu = Bf. 

The conjugate gradient method for SPD problems has several nice properties (good 
convergence rate, efficient three-term recursion, and minimization of the A-norm of the 
error at each step), some of which must be given up in order to generalize the method to 
nonsymmetric problems. For example, while GMRES attempts to maintain a minimiza- 
tion property and a good convergence rate, the three-term recursion must be sacrificed. 
Conjugate gradients on the normal equations maintains a minimization property as well 
as the efficient three-term recursion, but sacrifices convergence speed (the effective con- 
dition number is the square of the original system). Methods such as CGS and Bi-CGstab 
sacrifice the minimization property, but maintain good convergence speed and the effi- 
cient three-term recursion. For these reasons, methods such as CGS and Bi-CGstab have 
become the methods of choice in many applications that give rise to nonsymmetric prob- 
lems. Bi-CGstab has been shown to be more attractive than CGS in many situations due 
to the more regular convergence behavior [25]. In addition, Bi-CGstab does not require 
the application of the adjoint of the preconditioning operator, which can be difficult to 
implement in the case of some Schwarz methods. 

In §6, we shall use the preconditioned Bi-CGstab algorithm to accelerate nonsymmet- 
ric Schwarz methods. In a sequence of numerical experiments, we shall compare the 
effectiveness of this approach with unaccelerated symmetric and nonsymmetric Schwarz 
methods, and with symmetric Schwarz methods accelerated with conjugate gradients. 

3. Multiplicative Schwarz methods 

We develop a preconditioning theory of product algorithms which establishes suffi- 
cient conditions for producing SPD preconditioners. This theory is used to establish 
sufficient SPD conditions for multiplicative DD and MG methods. 
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3.1. A product operator. Consider a product operator of the form: 

E = I - BA = (I - B X A)(I - B A)(I - B X A) , 



(3.1) 



where B x , B and Bi are linear operators on %, and where A is, as before, an SPD 
operator on %. We are interested in conditions for B x , B and B x , which guarantee that 
the implicitly defined operator B is self-adjoint and positive definite and, hence, can be 
accelerated by using the conjugate gradient method. 

Lemma 3.1. Sufficient conditions for symmetry and positivity of operator B, implicitly 
defined by (3.1 ), are: 



(4) B non-negative on T-t . 

Proof. By Lemma 2.5, in order to prove symmetry of B, it is sufficient to prove that E 
is A-self-adjoint. By using (2.1), we get 



which follows from conditions 1 and 2. 

Next, we prove that (Bu,u) > 0, Vu G "H, u ^ 0. Since A is non-singular, this is 
equivalent to proving that (BAu, Au) > 0. Using condition 1, we have that 

(BAu,Au) = ({I-E)u,Au) 

= (u,Au) - ((I — BfA)(I — B A)(I — B x A)u, Au) 

= (u, Au) - ((/ - B A)(I - B x A)u, A(I - B x A)u) 

= (u, Au) - {(I - B x A)u, A(I - B x A)u) + (B w, w), 

where w = A(I — B x A)u. By condition 4, we have that (B w,w) > 0. Condition 3 
implies that ((/ — B x A)u, A(I — B x A)u) < (u, Au) for u ^ 0. Thus, the first two terms 
in the sum above are together positive, while the third one is non-negative, so that B is 
positive. □ 

Corollary 3.2. If B x = B\, then condition 3 in Lemma 3.1 is equivalent to p{I — B x A) < 
1. 

Proof. This follows directly from Lemma 2.1 and Lemma 2.5. □ 

3.2. Multiplicative domain decomposition. Given the finite-dimensional Hilbert space 
H, consider J spaces Hk, k — 1, . . . , J, together with linear operators I k 6 L(T-L k ,H), 
nuil(ifc) = {0}, such that IkHk Q H — Ylk=i ^krik- We also assume the existence 
of another space Ho, an associated operator I such that Jo^o Q an d some linear 
operators I k 6 L(H, Hk), k — 0, . . . , J. For notational convenience, we shall denote the 
inner-products on H k by (-, •) (without explicit reference to the particular space). Note 
that the inner-products on different spaces need not be related. 

In a domain decomposition context, the spaces Tik, k = 1, . . . , J, are typically as- 
sociated with local subdomains of the original domain on which the partial differential 
equation is defined. The space "H is then a space associated with some global coarse 
mesh. The operators Ik, k = 1, . . . , J, are usually inclusion operators, while I is an 



(1) Bi = Bf ; 

(2) B = 5 T ; 

(3) ||/-BiA|U<l; 



E* 



A- l E T A 

A' 1 (I - ABl){I - ABl){I - AB"()A 
(I-BfA)(I-B^A)(I-BfA) 
(I - B X A)(I - B A)(I - B X A) = E, 
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interpolation or prolongation operator (as in a two-level MG method). The operators 
I k ,k = 1, . . . , J, are usually orthogonal projection operators, while 7° is a restriction 
operator (again, as in a two-level MG method). 

The error propagator of a multiplicative DD method on the space % employing the 
subspaces I k T-L k has the general form [9]: 

E = I-BA=(I- IjRjI j A) •••(/- I RoI°A) •••(/- IjRjI j A) , (3.2) 

where R k and R k , k = 1, . . . , J, are linear operators on T-L k , and R is a linear operator 
on H . Usually the operators R k and R k are constructed so that R k ps A^ 1 and i4 

1, 



A fc \ where A fc is the operator defining the subdomain problem in T-L k . Similarly, R is 



constructed so that R ps A 1 . Actually, quite often R is a "direct solve", i.e., R = 
Aq 1 . The subdomain problem operator A k is related to the restriction of A to 7i k . We 
say that A k satisfies the Galerkin conditions or, in a finite element setting, that it is 
variationally defined when 

A k = I k AI k , I k = l T k . (3.3) 

Recall that the superscript "T" is to be interpreted as the adjoint in the sense of (2.2), i.e., 
with respect to the inner-products in H and T-L k . 

In the case of finite element, finite volume, or finite difference discretization of an 
elliptic problem, conditions (3.3) can be shown to hold naturally for both the matrices 
and the abstract weak form operators for all subdomains k = 1, . . . , J. For the coarse 
space T-L , often (3.3) must be imposed algebraically. 

Propagator (3.2) can be thought of as the product operator (3.1), by choosing 

i j 
I-B 1 A = Y[(I- I k R k I k A) , B = I RoI° , I — B X A = J](J - I k R k I k A) , 

k=J k=l 

where B\ and B\ are known only implicitly. (Note that we take the convention that the 
first term in the product appears on the left.) This identification allows for the use of 
Lemma 3.1 to establish sufficient conditions on the subdomain operators R k , R k and R 
to guarantee that multiplicative domain decomposition yields an SPD operator B. 

Theorem 3.3. Sufficient conditions for symmetry and positivity of the multiplicative do- 
main decomposition operator B, implicitly defined by (3.2), are: 

(1) I k = c k l T k , c k > , k = 0, • • • , J ; 

(2) R k = R T k , k = 1, • • • , J ; 

(3) Rq = Rq ; 



(4) U J k =i(I-hRkI k A) 

(5) Rq non-negative on T-L 



< 1 

A 



Proof. We show that the sufficient conditions of Lemma 3.1 are satisfied. First, we prove 
that Bi = Bf, which, by Lemma 2.5, is equivalent to proving that (J — B^A)* = 
(I — BiA). By using (2.1), we have 

\{{I-I k R k I k A)\ =A- 1 \T{{I-I k R k I k A)\ A = H(I-(I k ) T RT(I k ) T A), 

K k=l J \k=l J k=J 

which equals (/ — B^A) under conditions 1 and 2 of the theorem. The symmetry of B 
follows immediately from conditions 1 and 3; indeed, 

5 T = (/ #o/°) T = (I ) T Ro(hf = (co^Roico 1 ! ) = I RoI° = B . 
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By condition 4 of the theorem, condition 3 of Lemma 3.1 holds trivially. The theorem 
follows by realizing that condition 4 of Lemma 3.1 is also satisfied, since, 

(Bqu, u) = (I R I°u,u) = (R I°u,l£u) = c \R I u,I°u) > , Vw E U . 

□ 

Remark 3.4. Note that one sweep through the subdomains, followed by a coarse prob- 
lem solve, followed by another sweep through the subdomains in reversed order, gives 
rise an error propagator of the form (3.2). Also, note that no conditions are imposed on 
the nature of the operators A k associated with each subdomain. In particular, the theo- 
rem does not require that the variational conditions are satisfied. While it is natural for 
condition (3.3) to hold between the fine space and the spaces associated with each subdo- 
main, these conditions are often difficult to enforce for the coarse problem. Violation of 
variational conditions can occur, for example, when complex coefficient discontinuities 
do not lie along element boundaries on the coarse mesh ( we present numerical results for 
such a problem in §<5). The theorem also does not require that the overall multiplicative 
DD method be convergent. 

Remark 3.5. The results of the theorem apply for abstract operators on general finite- 
dimensional Hilbert spaces with arbitrary inner-products. They hold in particular for 
matrix operators on W\ equipped with the Euclidean inner-product, or the discrete L 2 
inner-product. In the former case, the superscript "T" corresponds to the standard ma- 
trix transpose. In the latter case, the matrix representation of the adjoint is a scalar 
multiple of the matrix transpose; the scalar may be different from unity when the adjoint 
involves two different spaces, and in the case of prolongation and restriction. This pos- 
sible constant in the case of the discrete L 2 inner-product is absorbed in the factor in 
condition 1. This allows for an easy verification of the conditions of the theorem in an 
actual implementation, where the operators are represented as matrices, and where the 
inner-products do not explicitly appear in the algorithm. 

Remark 3.6. Condition 1 of the theorem ( with — 1) for k = 1 , . . . , J is usually 
satisfied trivially for domain decomposition methods. For k = 0, it may have to be 
imposed explicitly. Condition 2 of the theorem allows for several alternatives which give 
rise to an SPD preconditioner, namely: (1) use of exact subdomain solvers (if Ak is a 
symmetric operator); (2) use of identical symmetric subdomain solvers in the forward 
and backward sweeps; (3) use of the adjoint of the subdomain solver on the second 
sweep. Condition 3 is satisfied when the coarse problem is symmetric and the solve 
is an exact one, which is usually the case. If not, the coarse problem solve has to be 
symmetric. Condition 4 in Theorem 3.3 is clearly a non-trivial one; it is essentially 
the assumption that the multiplicative DD method without a coarse space is convergent. 
Convergence theories for DD methods can be quite technical and depend on such things 
as the discretization, the subdomain number, shape, and size, and the regularity of the 
solution [5, 9, 27]. However, since variational conditions hold naturally between the 
fine space and each subdomain space for nearly any formulation of a DD method, very 
general convergence theorems can be derived, if one is not concerned about the actual 
rate of convergence. Using the Schwarz theory framework in any of [5, 9, 27], it can 
be shown that Condition 4 in Theorem 3.3 (convergence of multiplicative DD without a 
coarse space) holds if the variational conditions (3.3) hold, and if the subdomain solvers 
Rk are SPD. A proof of this result may be found for example in [14]. Condition 5 is 
satisfied for example when the coarse problem is SPD and the solve is exact. 
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Consider now the case when the subspaces together do not span the entire space, ex- 
cept when the coarse space is included. The above theorem can be applied with R = 0, 
and by viewing the coarse space as simply one of the spaces H k , k ^ 0. In this case, the 
error propagation operator E takes the form: 

I-BA = (I - IjRjI j A) ■ ■ ■ (I - hRJ 1 A)(I - hRJ 1 A) ■ ■ ■ (I - IjRjI j A) . (3.4) 

This leads to the following corollary. 

Corollary 3.7. Sufficient conditions for symmetry and positivity of the multiplicative 
domain decomposition operator B, implicitly defined by (3.4), are: 



(1) I k -- 

(2) R k 

(3) 



z R l 



k > 



C k > 

k = I,- 



hR k I k A) 



k = I,-- 
, J ; 
< 1 



J; 



Remark 3.8. Condition 3 is equivalent to requiring convergence of the overall multi- 
plicative Schwarz method. This follows from the relationship 

\\E\\ A =\\E*E\\ A = \\E\\ 2 A <1, 
where E = Y{ J k=1 (I - I k R k I k A). 

Remark 3.9. If in addition to conditions of the corollary, it holds that R\ = (7 1 AZi) -1 , 
i.e., it corresponds to an exact solve with a variationally defined subspace problem op- 
erator in the sense of (3.3), then 

(I - hR x I x A){I - hR^A) = 1- I X RJ X A, 

since I — I^I 1 AZi) -1 I\A is a projector. Therefore, space Hi (for example, the coarse 
space) needs to be visited only once in the application of (3.4). 

3.3. Multiplicative multigrid. Consider the Hilbert space "H, J spaces "H& together 
with linear operators Ik G L("H fc , H), null(J fc ) = 0, such that the spaces IkHk are nested 
and satisfy I {Hi C J 2 %2 Q • • • Q Ij-{Hj-i Q Hj = %• As before we denote the Th- 
inner-products by (•,•), since it will be clear from the arguments which inner-product is 
intended. Again, the inner-products are not necessarily related in any way. We assume 
also the existence of operators I k G L("H, Hk)- 

In a multigrid context, the spaces f-Lk are typically associated with a nested hierarchy 
of successively refined meshes, with Hi being the coarsest mesh, and Hj being the fine 
mesh on which the PDE solution is desired. The linear operators Ik are prolongation 
operators, constructed from given interpolation or prolongation operators that operate 
between subspaces, i.e., G L(Hk-i,Hk)- The operator I k is then constructed (only 
as a theoretical tool) as a composite operator 

h = rlii J jZ l 2 ---i k kXli k k + \ k = i,...,j-i. (3.5) 

The composite restriction operators I k , k — 1, . . . , J — 1, are constructed similarly from 
some given restriction operators I^ 1 G L(Hk, Hk-i)- 

The coarse problem operators A k are related to the restriction of A to H k - As in the 
case of DD methods, we say that A k is variationally defined or satisfies the Galerkin 
conditions when conditions (3.3) hold. It is not difficult to see that conditions (3.3) are 
equivalent to the following recursively defined variational conditions: 

A k = I k k+ iA k+l I k k + \ I k k+ i = (4 fc+1 ) T , (3.6) 
when the composite operators I k appearing in (3.3) are defined as in (3.5). 
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In a finite element setting, conditions (3.6) can be shown to hold in ideal situations, for 
both the stiffness matrices and the abstract weak form operators, for a nested sequence 
of successively refined finite element meshes. In the finite difference or finite volume 
method setting, conditions (3.6) must often be imposed algebraically, in a recursive fash- 
ion. 

The error propagator of a multiplicative V-cycle MG method is defined implicitly: 

E = I - BA = I - DjAj, (3.7) 

where Aj = A, and where operators D k , k = 2, . . . , J are defined recursively, 

/ - D k A k = (I- R k A k )(I - Jt 1 D fc _i/f-%)(/ - R k A k ), k = 2, . . . , J,(3.8) 
D x = R x . (3.9) 

Operators R k and R k are linear operators on T-L k , usually called smoothers. The linear 
operators A k E L(T-L k , H k ) define the coarse problems. They often satisfy the variational 
condition (3.6). 

The error propagator (3.7) can be thought of as an operator of the form (3.1) with 

B\ = Rj , B = Ij_iDj-\Ij 1 , Bi — Rj. 

Such an identification with the product method allows for use of the result in Lemma 3.1. 
The following theorem establishes sufficient conditions for the subspace operators R k , 
R k and A k in order to generate an (implicitly defined) SPD operator B that can be accel- 
erated with conjugate gradients. 

Theorem 3.10. Sufficient conditions for symmetry and positivity of the multiplicative 
multigrid operator B, implicitly defined by (3.7), (3.8), and (3.9), are 

(1) A k is SPD onH k , k = 2,...,J; 

(2) lt 1 = c k {lLi) T , c *> > k = 2,...,J; 

(3) R k = Rl, k = 2,...,J ; 

(4) R x = Rj ; 

(5) \\I-RjA\\ A <l,; 

(6) \\I-R k A k \\ Ak < 1, k = 2,...,J-l; 

(7) Ri non-negative on Hi . 

Proof. Since Rj = Rj, we have that Bi = Bf, which gives condition 1 of Lemma 3.1. 
Now, Bq is symmetric if and only if 

Bo = Ij ,/!/ ,/;/ 1 = {cfl-rfD^cjlUf = bi 

which holds under condition 2 and a symmetry requirement for We will prove 

that -Dj_i = D T j_ x by induction. First, Di = Df since Ri = Rj. By Lemma 2.5 and 
condition 1, D k is symmetric if and only if E k — I — D k A k is A fc -self-adjoint. By using 
(2.1), we have that 

El = A' 1 ((I - R k A k )(I - itiDk-iI^A^I - R k A k )f A k 

= A~\I - AlRl)(I - Al(lt-rDl_i(I k k _if)(I - A T k Rl)A k 

= (I- R T k A k )A k \l - Al{I^) T DU{lti) T )Ml - RU k ) 
= (I- R k A k )(I - (cJtJDT^I^A^I - R k A k ) , 

where we have used conditions 1, 2 and 3. Therefore, E% = E k , if D k ^i = D^_ 1 . Hence, 
the result follows by induction on k. 

Condition 3 of Lemma 3.1 follows trivially by condition 5 of the theorem. 
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It remains to verify condition 4 of Lemma 3.1, namely that B is non-negative. This 
is equivalent to showing that Dj_i is non-negative on Hj-x. This will follow again 
from an induction argument. First, note that D\ = R\ is non-negative on "Hi. Next, we 
prove that (D k v kl v k ) > 0, Vv k G T-i k , or, equivalently, since A k is non-singular, that 

(D k A k v k , A k v k ) > 0. So, for all v k G H k , 

(D k A k v k ,A k v k ) = (A k v k ,v k ) - (A k E k v k ,v k ) 
= (A k v k , v k ) 

-(A k (I - R k A k )(I - / fe fc _ 1 D fc _ 1 / fc fe - 1 A fe )(/ - R k A k )v k ,v k ) 
= (A k v k ,v k ) 

-(A k (I - I^D^J^A^I - R k A k )v k , (I - R k A k )v k ) 
= (A k v k , v k ) - (A k (I - R k A k )v k , (I - R k A k )v k ) 

+ (Akl^Dk^I*- 1 A k (I - R k A k )v k , (I - R k A k )v k ) 

= (V k , V k )A k ~ (S k V k , S k V k ) Ak + C^ 1 (D k _ 1 V k ^ 1 , V k ^i) 

where S k = I—R k A k and v k ^i = A k (I — R k A k )v k G H k -x. By condition 6, the first 
two terms in the above sum add up to a non-negative value. Hence, D k is non-negative 
if -Dfc_i is non-negative. Condition 4 of Lemma 3.1 follows. □ 

Corollary 3.11. If the fine grid smoother is symmetric, i.e., Rj = -Rj, then condition 5 
in Theorem 3.10 is equivalent to p{l — RjA) < 1. 

Proof. This follows directly from Corollary 3.2. □ 

Remark 3.12. The coarse grid operators A k , k = 2, . . . , J — 1, need only be SPD. They 
need not satisfy the Galerkin conditions (3.6). 

Remark 3.13. As noted earlier in Remark 3.5, the conditions and conclusions of the 
theorem can be interpreted completely in terms of the usual matrix representations of the 
multigrid operators. 

Remark 3.14. Condition 1 of the theorem requires that the coarse grid operators ( except 
for the coarsest one) be SPD. This is easily satisfied when they are constructed either by 
discretization or by explicitly using the Galerkin or variational condition. Condition 2 
requires restriction and prolongation to be adjoints, possibly multiplied by an arbitrary 
constant. Condition 3 of the theorem is satisfied when the number of pre-smoothing 
steps equals the number of post-smoothing steps, and in addition one of the following 
is imposed: (1) use of the same symmetric smoother for both pre- and post- smoothing; 
(2) use of the adjoint of the pre-smoothing operator as the post-smoother. Condition 4 
requires a symmetric coarsest mesh solver. When the coarsest mesh problem is SPD, 
the symmetry of Ri is satisfied when it corresponds to an exact solve ( as is typical for 
MG methods). Condition 5 is a convergence requirement on the fine space smoother. 
Condition 6 requires the coarse grid smoothers to be non-divergent. The nonnegativity 
requirement for Ri is a non-trivial one; however, if Ai is SPD, it is immediately satisfied 
when the operator corresponds to an exact solve. 

Theorem 3.10 applies to standard multigrid methods only. The conditions of the the- 
orem, and condition 5 in particular, cannot be satisfied in the cases of hierarchical basis 
multigrid methods [3], and multigrid methods with local smoothing on locally refined 
regions. The latter methods are covered in the following theorem, where the conditions 
that guarantee positivity of the preconditioner (conditions 5,6 and 7 in Theorem 3.10), 
are replaced by a convergence condition on the underlying iterative method. 
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Theorem 3.15. Sufficient conditions for symmetry and positivity of the multiplicative 
multigrid operator B, implicitly defined by (3.7), (3.8), and (3.9), are 

(1) A k is SPD onH k , k = 2, . . . , J ; 

(2) J*" 1 = c fc (i£_i) T , c k >0, k = 2, . . . , J ; 

(3) R k = RT, k = 2,...,J; 

(4) R, = BS[ ; 

(5) \\I-BA\\ A < 1. 

Proof. Positivity of B is proven easily by a contradiction argument. Symmetry follows 
from the proof of Theorem 3. 10. □ 

Requiring convergence of the underlying multigrid method is very restrictive; it is not 
a necessary condition. Positivity of B is satisfied if Aj(I — B A) < 1; no limit needs to 
be set on the magnitude of the negative eigenvalues. The above eigenvalue condition, 
however, does not seem to lead to conditions that are easily checked in practice. 

Remark 3.16. If variational conditions are satisfied on all levels, then there is a simple 
proof which shows that in addition to defining an SPD operator B, the conditions of 
Theorem 3.10 are sufficient to prove the convergence of the MG method itself. The result 
is as follows. 

Theorem 3.17. If in addition to the conditions for Theorem 3.10, it holds that A k = 
I k AI k , I k = I k , and Ri = A± , then the MG error propagator satisfies: 

p{E) < \\E\\ A < 1. 

Proof. Under the conditions of the theorem, the MG error propagator can be written 
explicitly as the product ([4, 19]): 

E=(I- IjR T jI T jA) •••(/- hRJ^A) •••(/- IjRjljA) . 

Since the coarse problem is solved exactly, and since variational conditions hold, the 
coarse product term is an A-orthogonal projector: 

/ - hRJ^A = 1- hihAlD'^lA = (I - hihAllyH^Af = (I - hRJ^Af. 

Therefore, we may define E — (I — IiRilfA) •••(/ — IjRjljA), and represent E as 
the product E = E*E. Now, since A is SPD, we have that: 

(AEv,v) = (AEv,Ev) > . 

Hence, E is A-non-negative. Under the conditions of the theorem, Lemma 3.1 implies 
that the preconditioner is SPD, and so by Lemma 2.9 it holds that \\E\\ A < 1. □ 

4. Additive Schwarz methods 

We now present an analysis of additive Schwarz methods. We establish sufficient 
conditions for additive algorithms to yield SPD preconditioner s. This theory is then 
employed to establish sufficient SPD conditions for additive DD and MG methods. 

4.1. A sum operator. Consider a sum operator of the following form: 

E = I — BA = I — uj(B + Bi)A, tu > , (4.1) 
where, as before, A is an SPD operator, and B and Bi are linear operators on H. 

Lemma 4.1. Sufficient conditions for symmetry and positivity of B, defined in (4.1), are 

(1) B x is SPD in H ; 

(2) Bq is symmetric and non-negative on % . 
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Proof. We have that B = cu(B + B{), which is symmetric by the symmetry of B and 
B\. Positivity follows since (B u, u) > and (BiU, u) > 0, Vu G H, u ^ 0. □ 

Remark 4.2. 77?e parameter u is usually required to make the additive method a con- 
vergent one. Its estimation is often nontrivial, and can be very costly. As was noted in 
Remark 2.13, the parameter u is not required when the linear additive method is used as 
a preconditioner in a conjugate gradients algorithm. This is exactly why additive multi- 
grid and domain decomposition methods are used almost exclusively as preconditioner s. 

4.2. Additive domain decomposition. As in §3.2, we consider the Hilbert space H, 
and J subspaces IkHk such that IkHk Q H — 2~2k=i IkHk- Again, we allow for the 
existence of a "coarse" subspace IqHq C %. 

The error propagator of an additive DD method on the space H employing the sub- 
spaces IkHk has the general form (see [27]): 

E = I - BA = I - oo(I R I° + hRJ 1 + ■■■ + IjRjI j )A. (4.2) 

The operators Rk are linear operators on Hk, constructed in such a way that Rk ~ A^ , 
where the Ak are the subdomain problem operators. Propagator (4.2) can be thought of 
as the sum method (4.1), by taking 

j 

Bo = IoRoI°, Bi = 2J hRk! k - 

k=l 

This identification allows for the use of Lemma 4.1 in order to establish conditions to 
guarantee that additive domain decomposition yields an SPD preconditioner. Before we 
state the main theorem, we need the following lemma, which characterizes the splitting 
of H into the subspaces IkHk in terms of a positive splitting constant Sq. 

Lemma 4.3. Given any v e H, there exists a splitting v = J2k=i Ik v k, Vk G Hk, and a 
constant So > 0, such that 

j 

52\\IkV k \\ 2 A<S \\v\\ 2 A . (4.3) 

fc=i 

Proof. Since J2k=i IkHk = H, we can construct subspaces 14 C Hk, such that 

j 

hVk n m = {0} , for k ^ I and H = I k V k . 

k=l 

Any v G H, can be decomposed uniquely as v — J2k=i^Vk, Vk G Vk- Define the 
projectors Q k G L(H, hVk) such that Q k v = IkVk- Then, 

j^\\hvkt A = Y^\\Qkv\\<jr \\Q k f A ht A . 

k=l k=l k=l 

Hence, the result follows with S = J2t=i HQklW- ^ 

Theorem 4.4. Sufficient conditions for symmetry and positivity of the additive domain 
decomposition operator B, defined in (4.2), are 

(1) I k = c k ll, c k >0,k = 0,...,J; 

(2) R k is SPD on H k , k = 1, . . . , J ; 

(3) Rq is symmetric and non-negative on Hq . 
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Proof. Symmetry of B and B\ follow trivially from the symmetry of R k and R , and 
from I k = c k I k . That B is non-negative on % follows immediately from the non- 
negativity of R on Ho- 

Finally, we prove positivity of B\. Define (only as a technical tool) the operators 
A k = I k AIk, k = 1, . . . , J. By condition 1, and the full rank nature of I k , we have that 
Ak is SPD. Now, since Rk is also SPD, the product RkAk is A^-SPD. Hence, there exists 
an tu > such that < tu < \i(R k A k ), k = 1, . . . , J. This is used together with (4.3) 
to bound the following sum, 

j j 
^Z c k\ R k lv k^ k ) = ^2c k 1 (A k A- 1 R- 1 Vk,Vk) 



k=i k=i 



\ -[, , . {A-kA k 1 R k v k , v k ) \ 

^ l^ c k {A k v k ,v k )max — r < > c k u {A k v k ,v k ) 

^ v k ^o (A k v k ,v k ) ^ 



y^o; i (Ai k v k , i k v k ) = y^uj ^i/fc^fciii < ( — ) \\v\\ A . 



A;=l fc=l 

with v = J2k=i Ik v k- We can now employ this result to establish positivity of Bi. First 
note that 



./ 



MB 



(A;,*) = ^(^,4^) = ^(4 T ^, Wfc ) = J2( R d /2 llAv^cl^Vk) ■ 
k=i fc=i fc=i 

By using the Cauchy-Schwarz inequality first in the /4-inner-product and then in IR J , 
we have that 

(3 \ V2 / j X 1/2 

IMft < ^~ V 7 V) ^(RkC^llAv, c^I T k Av) 



< 



,k=l / \k=l 

1/2 

'0 



' o \ 1/2 
On 



^0 



OQ 



Finally, division by ||f ||a an d squaring yields 

U II. .112 



D 



□ 



Remark 4.5. Condition 1 is naturally satisfied for k = 1, . . . , J, with c k = 1, since the 
associated I k and I k are usually inclusion and orthogonal projection operators (which 
are natural adjoints when the inner-products are inherited from the parent space, as 
in domain decomposition). The fact that 1° = CqIq needs to be satisfied explicitly. 
Condition 2 requires the use of SPD subdomain solvers. The condition will hold, for 
example, when the subdomain solve is exact and the subdomain problem solver is SPD. 
(The latter is naturally satisfied by condition 1 and the full rank nature of I k .) Finally, 
condition 3 is nontrivial, and needs to be checked explicitly. The condition holds when 
the coarse space problem operator is SPD and the solve is exact. Note that variational 
conditions are not needed for the coarse space problem operator. 
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Consider again the case when the subspaces together do not span the entire space, 
except when the coarse space is included. The above theorem applies immediately with 
Ro = 0, where now the coarse space is taken to be any one of the spaces k ^ 0. The 
error propagator takes the form 

I-BA = I- uj{hRJ l + I 2 R 2 I 2 + ■■■ + IjRjI j )A. (4.4) 

This leads to the following corollary. 

Corollary 4.6. Sufficient conditions for symmetry and positivity of the additive domain 
decomposition operator B, defined in (4.4), are 

(1) I k = c k I%, c k >Q,k = l,...,J; 

(2) R k is SPD on 7-L k , k = 1, . . . , J . 

A3. Additive multigrid. As in §3.3, given are the Hilbert space H, and J — 1 nested 
subspaces I k ?i k such that I {Hi C 1{H 2 Q • • • Q Ij-iHj-i Q Hj = H . The operators 
I k and I k are the usual linear operators between the different spaces, as in the previous 
sections. 

The error propagator of an additive MG method is defined explicitly: 

E = I - BA = I - ooihRJ 1 + I 2 R 2 I 2 + ■■■ + Ij-iRj^J 1 - 1 + Rj)A. (4.5) 

This can be thought of as the sum method analyzed earlier, by taking 

J -1 

Bo = IkRkI k , Bi = Rj . 

k=l 

This identification allows for the use of Lemma 4.1 to establish sufficient conditions 
to guarantee that additive MG yields an SPD preconditioned 

Theorem 4.7. Sufficient conditions for symmetry and positivity of the additive multigrid 
operator B, defined in (4.5), are: 

(1) I k = c k ll , c k > , k = 1, . . . , J - 1 ; 

(2) Rj is SPD in H; 

(3) R k is symmetric non-negative in % k , k — 1, . . . , J — 1. 

Proof. Symmetry of B and B\ is obvious. B% is positive by condition 2. Non-negativity 
of Bo follows from 

j-i j-i 
(B u, u) = ^2(I k R k (c k I k ) T u, u) = ^ c k (R k I k u, ilu) > 0, \fuen,u^0. 

k=l k=l 

□ 

Remark 4.8. Condition 1 of the theorem has to be imposed explicitly. Conditions 2 and 3 
require the smoothers to be symmetric. The positivity of Rj is satisfied when the fine grid 
smoother is convergent, although this is not a necessary condition. The non-negativity 
of R k , k < J, has to be checked explicitly. When the coarse problem operators are SPD, 
this condition is satisfied, for example, when the smoothers are non-divergent. Note that 
variational conditions for the subspace problem operators are not required. 

Theorem 4.7 is applicable to the standard multigrid case, i.e., the case where the fine 
mesh smoother operates on the entire fine mesh. As in §3.3, a different set of conditions 
is to be derived to cover the cases of (additive) hierarchical basis preconditioners [28], 
and additive multilevel methods with smoothing in local refinement regions only. The 
latter cases are treated most easily by loosening the restriction of nestedness of the spaces 
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Ik'Hk- With "H fe interpreted as the domain space of the smoother R k , the theory becomes 
identical to that of the additive domain decomposition case. Sufficient conditions for the 
additive operator B to be SPD are then similar to the conditions of the additive domain 
decomposition method in Corollary 4.6 

5. TO SYMMETRIZE OR NOT TO SYMMETRIZE 

The following lemma illustrates why symmetrizing is a bad idea for linear methods. It 
exposes the convergence rate penalty incurred by symmetrization of a linear method. 

Lemma 5.1. For any E 6 L(H, H), it holds that: 

p(EE) < \\EE\\ A < \\E\\ 2 A = \\EE*\\ A = p{EE*). 

Proof. The first and second inequalities hold for any norm. The first equality follows 
from Lemma 2.8, and the second follows from Lemma 2.1. □ 

Note that this is an inequality not only for the spectral radii, which is only an asymp- 
totical measure of convergence, but also for the A-norms of the nonsymmetric and sym- 
metrized error propagators. The lemma illustrates that one may actually see the differing 
convergence rates early in the iteration as well. 

Based on this lemma, and Corollary 2.16, we conjecture that when symmetrization 
of a linear method is required for its use as a preconditioner, the best results will be 
obtained by enforcing only a minimal amount of symmetry. This conjecture can be made 
more clear by considering the two linear methods with error propagators E\ = EE 
and E 2 = EE*, and the effectiveness as preconditioners of the symmetrized product 
operators E\E{ and E 2 E\. Using Lemma 2.1 and Lemma 5.1, we find immediately that 

ll-E'i-E'i IU = H-Eilll < II-^Ha = II-^-EjIU- (5.1) 

While both operators are symmetric, the operator E\E* has been symmetrized "mini- 
mally" in the sense that the individual terms E x making up the product are themselves 
nonsymmetric. On the other hand, both of the terms E 2 making up the operator E 2 E\ 
are completely symmetric. 

As a result of inequality (5.1) and Corollary 2.16, the bound for the condition number 
of the preconditioner associated with E\E* is less than the corresponding bound for the 
preconditioner associated with E 2 E\. Hence, the "less-symmetric" E±E* would likely 
produce a better preconditioner than the "more- symmetric" E 2 E^. 

6. Numerical results 

We present numerical results obtained by using multiplicative and additive finite- 
element-based DD and MG methods applied to two test problems, and we illustrate the 
theory of the preceding sections. 

6.1. Example 1. Violation of variational conditions can occur in DD and MG methods 
when, for example, complex coefficient discontinuities do not lie along element bound- 
aries on coarse meshes. An example of this occurs with the following test problem. The 
Poisson-Boltzmann equation describes the electrostatic potential of a biomolecule lying 
in an ionic solvent (see, e.g., [7] for an overview). This nonlinear elliptic equation for 
the dimensionless electrostatic potential u(r) has the form: 

-V • (e(r)Vu(r)) + R 2 sinh(u(r)) = ( J ^ z l 5{r - n), r e M 3 , u(oo) = . 
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Table 1 . Normalized operation counts per iteration, Example 1. 



Method 


UNACCEL 


CG 


Bi-CGstab 


multiplicative MG 


1.0 


1.4 


2.6 


additive MG 


.95 


1.3 


2.5 


multiplicative DD 


3.5 


3.8 


7.5 


additive DD 


3.1 


3.4 


6.7 



The coefficients appearing in the equation are discontinuous by several orders of magni- 
tude. The placement and magnitude of atomic charges are represented by source terms 
involving delta-functions. Analytical techniques are used to obtain boundary conditions 
on a finite domain boundary. 

We will compare several MG and DD methods for a two-dimensional, linearized 
Poisson-Boltzmann problem, modeling a molecule with three point charges. The sur- 
face of the molecule is such that the discontinuities do not align with the coarsest mesh 
or with the subdomain boundaries. Beginning with the coarse mesh shown on the left 
in Figure 1, we uniformly refine the initial mesh of 10 elements (9 nodes) four times, 
leading to a fine mesh of 2560 elements (1329 nodes). Piecewise linear finite elements, 
combined with one-point Gaussian quadrature, are used to discretize the problem. The 
three coarsest meshes used to formulate the MG methods are given in Figure 1 . For the 
DD methods, the subdomains, corresponding to the initial coarse triangulation, are given 
a small overlap of one fine mesh triangle. The DD methods also employ a coarse space 
constructed from the initial triangulation. Figure 2 shows three overlapping subdomains 
overlaying the initial coarse mesh. 

Computed results are presented in Tables 2 to 5. Given for each experiment is the 
number of iterations required to satisfy the error criterion (reduction of the A-norm of 
the error by 10~ 10 ). We report results for the unaccelerated, CG-accelerated, and Bi- 
CGstab-accelerated methods. Since the cost of one iteration differs for each method, 
Table 1 gives the operation counts per iteration, normalized by the cost of a single multi- 
grid iteration. For the MG operation counts, two smoothing iterations by lexicographic 
Gauss-Seidel (one pre- and one post-smoothing) are used. The DD operation counts are 
for methods employing two sweeps through the subdomains, each approximate subdo- 
main solve consisting of four sweeps of a Gauss-Seidel iteration. 

Table 1 shows that multiplicative MG is slightly more costly than additive MG, since 
additive MG requires the computation of the residual only on the finest level. Similarly, 
multiplicative DD is somewhat more costly than additive DD, due to the need to up- 
date boundary information (recompute the residual) after the solution of each subdomain 
problem. Table 1 should not be used to compare MG and DD methods for efficiency. 
Similar experiments [15] with more carefully optimized DD and MG methods show DD 
to be often competitive with MG for difficult elliptic equations such as those with discon- 
tinuous coefficients, although there may be some debate as to which approach is more 
effective on parallel computers [24] . 

Multiplicative multigrid. The results for multiplicative V-cycle MG are presented in Ta- 
ble 2. Each row corresponds to a different smoothing strategy, and is annotated by 
(uijVz), with V\. pre-smoothing strategy, and v 2 : post-smoothing strategy. An "f" in- 
dicates the use of a single forward Gauss-Seidel sweep, while a "b" denotes the use of 
the adjoint of the latter, i.e., a backward Gauss-Seidel sweep. u 2 ) = (//, fb), for ex- 
ample, corresponds to two forward Gauss-Seidel pre-smoothing steps, and a symmetric 
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Figure 1 . Example 1: Nested finite element meshes for MG. 




Figure 2. Example 1: Overlapping subdomains for DD. 



(forward/backward) post- smoothing step. Two series of results are given. For the first 
set, we explicitly imposed the Galerkin conditions when constructing the coarse opera- 
tors. In this case, the multigrid algorithm is guaranteed to converge by Theorem 3.17. 
In the second series of tests, corresponding to the numbers in parentheses, the coarse 
mesh operators are constructed using standard finite element discretization. In that case, 
Galerkin conditions are not satisfied everywhere due to coefficient discontinuities ap- 
pearing within coarse elements; hence, the MG method may diverge (DIV). 

The unaccelerated MG results clearly illustrate the symmetry penalty discussed in 
§5. The nonsymmetric methods are always superior to the symmetric ones (the cases 
(f,b), (ff,bb), and (fb,fb)). Note that minimal symmetry (ff,bb) leads to a better conver- 
gence than maximal symmetry (fb,fb). The correctness of Lemma 5.1 is illustrated by 
noting that two iterations of the (f,0) strategy are actually faster than one iteration of 
the (f,b) strategy; also, compare the (ff,0) strategy to the (ff,bb) one. CG-acceleration 
leads to a guaranteed reduction in iteration count for the symmetric preconditioners (see 
Lemma 2.17). We observe that the unaccelerated method need not be convergent for 
CG to be effective (recall Remarks 2.13 and 4.2, and the (f,b) result). CG appears to 
accelerate also some non-symmetric linear methods. Yet, it seems difficult to predict 
failure or success beforehand in such cases. The most robust method appears to be the 
Bi-CGstab method. The number of iterations with this method depends only marginally 
on the symmetric or nonsymmetric nature of the linear method. Note the tendency to 
favor the nonsymmetric V-cycle strategies. Overall, the fastest method proves to be the 
Bi-CGstab-acceleration of a (very nonsymmetric) V(l,0)-cycle. 

Multiplicative domain decomposition. Some numerical results for multiplicative DD 
with different subdomain solvers, and different subdomain sweeps are given in Table 3. 
In the column "forw", the iteration counts reported were obtained with a single sweep 
though the subdomains on each multiplicative DD iteration. The other columns corre- 
spond to a symmetric forward/backward sweep or to two forward sweeps. Four different 
subdomain solvers are used: an exact solve, a symmetric method consisting of two sym- 
metric Gauss-Seidel iterations, a nonsymmetric method consisting of four Gauss-Seidel 
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iterations, and, finally, a method using four forward Gauss-Seidel iterations in the for- 
ward subdomain sweep and using their adjoint, i.e., four backward Gauss-Seidel itera- 
tions, in the backward subdomain sweep. The latter leads to an symmetric iteration; see 
Remark 3.4. Note that the cost of the three inexact subdomain solvers is identical. 

Although apparently not as sensitive to operator symmetries as MG, the same con- 
clusions can be drawn for DD as for MG. In particular, the symmetry penalty is seen 
for the pure DD results. Lemma 5.1 is confirmed since two iterations in the column 
"forw" are always more efficient than one iteration of the corresponding symmetrized 
method in column "forw/back". The CG results indicate that using minimal symmetry 
(the "adjointed" column) is a more effective approach than the fully symmetric one (the 
"symmetric" column). Again, the most robust acceleration is the Bi-CGstab one. 

Additive multigrid. Results obtained with an additive multigrid method are reported in 
Table 4. The number and nature of the smoothing strategy is given in the first column of 
the table. 

In the case of an unaccelerated additive method, the selection of a good damping pa- 
rameter is crucial for convergence of the method. We did not search extensively for an 
optimal parameter; a selection of to = 0.45 seemed to provide good results in the case 
when the coarse problem is variationally defined. No w-value leading to satisfactory con- 
vergence was found in the case when the coarse problem is obtained by discretization. In 
the case of CG acceleration the observed convergence behavior was completely indepen- 
dent of the choice of u; see Remark 3.4. The symmetric methods (y = fb, ffbb, fbfb) 
are accelerated very well. Some of the nonsymmetric methods are accelerated too, es- 
pecially when the number of smoothing steps is sufficiently large. In the case of Bi- 
CGstab-acceleration, there appeared to be a dependence of convergence on u (only with 
use of non-variational coarse problem). In that case we took to — 1. The overall best 
method appears to be the Bi-CGstab acceleration of the nonsymmetric multigrid method 
with a single forward Gauss-Seidel sweep on each grid-level. 

Additive domain decomposition. The results for additive DD are given in Table 5. The 
subdomain solver is either an exact solver, a symmetric solver based on two symmet- 
ric (forward/backward) Gauss-Seidel sweeps, or a nonsymmetric solver based on four 
forward Gauss-Seidel iterations. 

No value of to was found that led to satisfactory convergence of the unaccelerated 
method. CG-acceleration performs well when the linear method is symmetric; it per- 
forms less well for the nonsymmetric method. Again, the best overall method is the 
Bi-CGstab-acceleration of the nonsymmetric additive solver. 

6.2. Example 2. The second test problem is the Laplace equation on a semi-adapted 
L-shaped domain, with Dirichlet boundary conditions chosen in such a way that the 
equation has the following solution (in polar coordinates): 

u(r,6) = Vr sin(0/2) , 

where the re-entrant corner in the domain is located at the origin. Note that the one-point 
Gaussian quadrature rule which we employ to construct the stiffness matrix entries is an 
exact integrator here. Hence, the variational conditions (3.3) hold automatically between 
the fine space and all subdomain and coarse spaces for both the MG and the DD methods. 

Figure 3 shows a nested sequence of uniform mesh refinements used to formulate the 
MG methods. (A total of 5 mesh levels is used in the computation.) Figure 4 shows sev- 
eral overlapping subdomains constructed from a piece of the fine mesh of 9216 elements 
(4705 nodes) overlaying the initial coarse mesh of 36 elements (25 nodes). 
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Table 2. Example 1: Multiplicative MG with variational (discretized) 
coarse problem 





UNACCEL 


CG 


Bi-CGstab 


f 


65 (DIV) 


>100 (>100) 


14 (16) 


f b 


55 (DIV) 


16 (18) 


10 (15) 


f f 


40 (31) 


30 (>100) 


9 (9) 


ff 


39 (48) 


>100 (»100) 


8 (10) 


fb 


53 (DIV) 


>100 (>100) 


10 (11) 


ff 


39 (29) 


29 (>100) 


8 (9) 


fb 


53 (DIV) 


17 (99) 


10 (12) 


fb fb 


34 (27) 


12 (13) 


8 (8) 


ff bb 


28 (18) 


11 (11) 


7 (7) 


ff ff 


24 (15) 


12 (12) 


6 (6) 


fff f 


24 (15) 


17 (27) 


6 (6) 


ffff 


25 (17) 


>100 (»100) 


7 (6) 



Table 3. Example 1: Multiplicative DD with variational (discretized) 
coarse problem 



Accel. 


subdomain solve 


forw 


forw/back 


forw/forw 


UNACCEL 


exact 


40 (42) 


38 (39) 


20 (21) 


symmetric 


279 (282) 


146 (149) 


140 (141) 


adjointed 




110 (112) 


102 (103) 


nonsymmetric 


189 (191) 


102 (104) 


95 (96) 


CG 


exact 


>500 (>500) 


13 (13) 


20 (20) 


symmetric 


140 (56) 


24 (24) 


29 (27) 


adjointed 




21 (21) 


25 (26) 


nonsymmetric 


135 (83) 


22 (23) 


28 (28) 


Bi-CGstab 


exact 


9 (9) 


9 (9) 


6 (6) 


symmetric 


23 (23) 


17 (16) 


16 (16) 


adjointed 




14 (14) 


14 (13) 


nonsymmetric 


19 (20) 


13 (13) 


13 (13) 



Table 4. Example 1: Additive MG with variational (discretized) coarse problem 



V 


UNACCEL 


CG 


Bi-CGstab 


f 


175 (>1000) 


>100 (>100) 


23 (52) 


ff 


110 (>1000) 


119 (168) 


19 (43) 


fb 


146 (>1000) 


34 (54) 


23 (49) 


ffff 


95 (>1000) 


28 (67) 


17 (37) 


ffbb 


100 (>1000) 


27 (47) 


17 (34) 


fbfb 


95 (>1000) 


28 (48) 


20 (43) 
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Table 5. Example 1: Additive DD with variational (discretized) coarse problem 



subdomain solve 


UNACCEL 


CG 


Bi-CGstab 


exact 


»1000 (>1000) 


34 (34) 


25 (27) 


symmetric 


»1000 (>1000) 


57 (57) 


50 (49) 


nonsymmetric 


>1000 (>1000) 


69 (65) 


38 (41) 




Figure 3. Example 2: Nested finite element meshes for MG. 




Figure 4. Example 2: Overlapping subdomains for DD. 

Multiplicative Methods. The results for multiplicative MG are given in Table 6, whereas 
the results for multiplicative DD are given in Table 7. The results are similar to those 
for Example 1; in particular, imposing minimal symmetry is the most effective CG- 
accelerated approach to the problem. Employing the least symmetric linear method alone 
is the most effective linear method, and the same nonsymmetric linear method yields the 
most effective Bi-CGstab-accelerated approach. 

Additive Methods. As for Example 1, in the case of the unaccelerated additive methods 
the selection of the damping parameter was crucial for convergence of the methods. We 
did not search extensively for an optimal parameter; a selection of to = 0.45 seemed 
to provide acceptable results for DD. Note that improved convergence behavior might 
be obtained by allowing different u values for each subdomain solver (this will not be 
further investigated here). No satisfactory value for u was found for additive MG. In the 
case of CG acceleration, the observed convergence behavior was completely independent 
of the choice of to. The results for additive MG are given in Table 8, whereas the results 
for additive DD are given in Table 9. The effect of the symmetry of the linear method's 
error propagator on its convergence, and on the convergence behavior of CG and Bi- 
CGstab, was as for Example 1 . 

7. Concluding remarks 

In this paper, we developed framework for establishing sufficient conditions which 
guarantee that abstract multiplicative and additive Schwarz algorithms to yield self- 
adjoint positive definite preconditioners. We then analyzed four specific methods: MG 
and DD methods, in both their additive and multiplicative forms. In all four cases, we 
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used the general theory to establish sufficient conditions that guarantee the resulting 
preconditioner is SPD. As discussed in Remarks 3.4, 3.12, 4.5, and 4.8, the sufficient 
conditions for the theory, in the case of all four methods, are easily satisfied for non- 
variational, and even non-convergent methods. The analysis shows that by simply taking 
some care in the way a Schwarz method is formulated, one can guarantee that the method 
is convergent when accelerated with the conjugate gradient method. 

We also investigated the role of symmetry in linear methods and preconditioners. A 
certain penalty lemma (Lemma 5.1) was stated and proved, illustrating why symmetriz- 
ing is actually a bad idea for linear methods. It was conjectured that enforcing minimal 
symmetry in a linear preconditioner achieves the best results when combined with the 
conjugate gradient method, and our numerical examples illustrate this behavior almost 
uniformly. A sequence of experiments with two non-trivial test problems showed that the 
most efficient approach may be to abandon symmetry in the preconditioner altogether, 
and to employ a nonsymmetric solver such as Bi-CGstab. While acceleration with CG 
was strongly dependent on the symmetric nature of the preconditioner, Bi-CGstab al- 
ways converged rapidly. In addition, BiCGstab appeared to benefit from the behavior 
predicted by Lemma 5.1, namely that a nonsymmetric linear preconditioner should have 
better convergence properties than its symmetrized form. 
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Table 6. Example 2: Multiplicative MG 



25 





V\ 


^2 


UNACCEL 


CG 


Bi-CGstab 






f 





33 


>100 


11 






f 


b 


23 


12 


7 






f 


f 


19 


22 


6 






ff 





21 


>100 


7 






fb 





25 


»100 


8 









ff 


20 


42 


7 









fb 


23 


17 


8 






fb 


fb 


16 


9 


6 






ff 


bb 


15 


9 


5 






ff 


ff 


14 


9 


5 






fff 


f 


14 


12 


5 






ffff 





16 


36 


5 




Table 7 


. Example 2: Multiplicative DD 




Accel. 


subdomain solve 


forw 


forw/back 


forw/forw 


UNACCEL 


exact 


73 


60 


37 




symmetric 


402 


205 


207 




adjointed 




153 


146 




nonsymmetric 


267 


144 


134 


CG 




exact 


116 


17 


17 




symmetric 


164 


37 


38 




adjointed 




32 


33 




nonsymmetric 


121 


31 


32 


Bi-CGstab 


exact 


11 


11 


7 




symmetric 


37 


25 


26 




adjointed 




22 


23 




nonsymmetric 


27 


21 


21 



Table 8. Example 2: Additive MG 



V 


UNACCEL 


CG 


Bi-CGstab 


f 


91 


>1000 


21 


ff 


62 


31 


16 


fb 


74 


29 


18 


ffff 


126 


25 


14 


ffbb 


136 


27 


15 


fbfb 


98 


27 


15 



Table 9. Example 2: Additive DD 



subdomain solve 


UNACCEL 


CG 


Bi-CGstab 


exact 


>1000 


42 


29 


symmetric 


»1000 


86 


56 


nonsymmetric 


>1000 


82 


49 



